#include "base.h"

int main(int argc,char **argv)
{

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@///

	  //////////////////
	  // Input   ///////
	  //////////////////

	  int Nr=2;
	  char kappa_sea[1024];
	  int conf_ini, conf_end, jump;
	  int Nconf, Nex;
	  int Nmass;
	  char dir_input[1024];
	  char dir_output[1024];

	

	  ////// Parsing command line

    	  for(int i=1; i<argc; i++) {
	    if(strncmp(argv[i],"-kappa_sea",10)==0)  sscanf(argv[i],"-kappa_sea=%s",kappa_sea);
	    if(strncmp(argv[i],"-conf_ini",9)==0)  sscanf(argv[i],"-conf_ini=%d",&conf_ini);
	    if(strncmp(argv[i],"-jump",5)==0)  sscanf(argv[i],"-jump=%d",&jump);
	    if(strncmp(argv[i],"-Nclust",7)==0)  sscanf(argv[i],"-Nclust=%d",&Nclust);
	    if(strncmp(argv[i],"-Nex",4)==0)  sscanf(argv[i],"-Nex=%d",&Nex);
	    if(strncmp(argv[i],"-Nmass",6)==0)  sscanf(argv[i],"-Nmass=%d",&Nmass);
            if(strncmp(argv[i],"-dir_input",10)==0)  sscanf(argv[i],"-dir_input=%s",dir_input);
            if(strncmp(argv[i],"-dir_output",11)==0)  sscanf(argv[i],"-dir_output=%s",dir_output);
	    if(strncmp(argv[i],"-T",2)==0)  sscanf(argv[i],"-T=%d",&T); 
            if(strncmp(argv[i],"-L",2)==0)  sscanf(argv[i],"-L=%d",&L); 

    	  }

	  Nconf=Nclust*Nex;
	  conf_end=conf_ini+jump*(Nconf-1);
          printf("==============================\n");
          printf("conf_ini=%d, conf_end=%d jump=%d\nNconf=%d, Nclust=%d, Nex=%d\n",conf_ini, conf_end,jump,Nconf,Nclust,Nex);
          printf("Nmass=%d\n",Nmass);
          printf("L=%d T=%d\n",L,T);
          printf("Nr=%d \n",Nr);
          printf("dir_input=%s \n",dir_input);
          printf("dir_output=%s \n",dir_output);
          printf("==============================");
	
   /////////////////////////
  //// Momenta set
  //////////////////////////////
  char momentafile[1024]; sprintf(momentafile,"../INPUT/mom.dat");
  read_momenta_data(momentafile);
  double *nn=(double*)malloc(sizeof(double)*(Nmom));
  compute_nn(nn);



	  //////////////////////////////////
          // IO files
          ///////////////////////////////

          char s0ff_base_input[1024];
          char s0ft_input[1024];
          char s0ff_base_output[1024];
          char s0ff_output[1024];
          char z_base_output[1024];
          char z_output[1024];
	  char Lambda_base_output[1024];
	  char Lambda_output[1024];

	  char file[1024];

	  char gammaprint[8];
        
          sprintf(s0ff_base_input,"%s/s0ft",dir_input);

          sprintf(z_base_output,"%s/Z_ALLP/Z",dir_output);
          sprintf(s0ff_base_output,"%s/FF_ALLP/ff",dir_output);
          sprintf(Lambda_base_output,"%s/LAMBDA_ALLP/Lambda",dir_output);


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//

           /////////////////////////////
	   // Read Fourier transform
	   //////////////////////////////
	
          //def
	  colorcolorspinspin ****s0FT_r=malloc(sizeof(colorcolorspinspin***)*(Nconf));
          for (int iconf=0; iconf<Nconf; iconf++) s0FT_r[iconf]=malloc(sizeof(colorcolorspinspin**)*(Nmass));
          for (int iconf=0; iconf<Nconf; iconf++)  for(int imass=0; imass<Nmass; imass++) s0FT_r[iconf][imass]=malloc(sizeof(colorcolorspinspin*)*(Nmom)); 
	  for (int iconf=0; iconf<Nconf; iconf++)  for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++) s0FT_r[iconf][imass][imom]=malloc(sizeof(colorcolorspinspin)*(Nr));
	
	  int conf_to_read;
	  FILE *fin=NULL;


           printf("====================\n");
           printf("READING...\n");

           for (int r=0; r<Nr; r++)
	   {
           for (int iconf=0; iconf<Nconf; iconf++)
 	   {
	           conf_to_read=conf_ini+iconf*jump;
                   //printf("  ... conf %d\n",conf_to_read);
                   sprintf(s0ft_input,"%s%d_%s_%03d.out",s0ff_base_input,r,kappa_sea,conf_to_read);
		   sprintf(s0ft_input,"%s%d_%04d.out",s0ff_base_input,r,conf_to_read);
                   fin=fopen(s0ft_input,"rb+");
                   if(fin == (FILE *)NULL) {
                        fprintf(stderr,"Unable to open file %s for reading\n",s0ft_input);
			sprintf(s0ft_input,"%s%d_%04d.out",s0ff_base_input,r,conf_to_read);
			fin=fopen(s0ft_input,"rb+");
                        if(fin == (FILE *)NULL) return EXIT_FAILURE;
			else printf("Instead opening %s \n",s0ft_input);
                   }
	           for(int imass=0; imass<Nmass; imass++)
                   {
                  // printf("im=%d\n",imass);
		   for (int imom=0; imom<Nmom; imom++)
			for (int ic_ext=0; ic_ext<3; ic_ext++) for (int is_ext=0; is_ext<4; is_ext++) for (int ic_int=0; ic_int<3; ic_int++) for (int is_int=0; is_int<4; is_int++)
			fread(&s0FT_r[iconf][imass][imom][r][ic_int][ic_ext][is_int][is_ext],sizeof(complex),1,fin);
		   }
		   fclose(fin);
    	
   	   }
	   }
	   printf("====================\n\n\n");

	    ////////////////////////////////////
            /////         CLUSTERIZE      //////
	    ////////////////////////////////////
	    printf("====================\n");
            printf("CLUSTERIZE\n");

	   complex   aux_conf[Nconf], aux_clust[Nclust+1];

          //def
          colorcolorspinspin ****s0FT=malloc(sizeof(colorcolorspinspin***)*(Nclust+1));
          for (int iclust=0; iclust<Nclust+1; iclust++) s0FT[iclust]=malloc(sizeof(colorcolorspinspin**)*(Nmass));
          for (int iclust=0; iclust<Nclust+1; iclust++)  for(int imass=0; imass<Nmass; imass++) s0FT[iclust][imass]=malloc(sizeof(colorcolorspinspin*)*(Nmom));
	  for (int iclust=0; iclust<Nclust+1; iclust++)  for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++) s0FT[iclust][imass][imom]=malloc(sizeof(colorcolorspinspin)*(Nr));

	  for (int r=0; r<Nr; r++)
	    for(int imass=0; imass<Nmass; imass++)  for (int imom=0; imom<Nmom; imom++)
            {
		for (int ic_ext=0; ic_ext<3; ic_ext++) for (int is_ext=0; is_ext<4; is_ext++) for (int ic_int=0; ic_int<3; ic_int++) for (int is_int=0; is_int<4; is_int++)
		{	
	               for(int iconf=0;iconf<Nconf;iconf++) aux_conf[iconf]=s0FT_r[iconf][imass][imom][r][ic_int][ic_ext][is_int][is_ext];
                       conf2clust(aux_conf, aux_clust,Nconf,Nclust,Nex);
                       for(int iclust=0;iclust<Nclust+1;iclust++) s0FT[iclust][imass][imom][r][ic_int][ic_ext][is_int][is_ext]=aux_clust[iclust];
		}
            }
	   
            printf("====================\n\n\n");




//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//

	/////////////////////////////////////////////////////
        // FORM FACTOR OF THE FOURIER-TRANSFORMED PROPAGATOR
	//////////////////////////////////////////////////////

        printf("====================\n");
        printf("FORM FACTOR\n");

	//def
	complex *****s0_ff=malloc(sizeof(complex****)*(Nclust+1));
        for (int iclust=0; iclust<Nclust+1; iclust++) s0_ff[iclust]=malloc(sizeof(complex***)*(Nmass));
        for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) s0_ff[iclust][imass]=malloc(sizeof(complex**)*(Nmom));
	 for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++) s0_ff[iclust][imass][imom]=malloc(sizeof(complex*)*(Nr));
	for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++) for (int r=0; r<Nr; r++) s0_ff[iclust][imass][imom][r]=malloc(sizeof(complex)*(3));
	

	 for (int imom=0; imom<Nmom; imom++)
	 {
		        for (int iclust=0; iclust<Nclust+1; iclust++) for (int imass=0; imass<Nmass; imass++) for (int r=0; r<Nr; r++)
			compute_s0tw_ff(s0FT[iclust][imass][imom][r],s0_ff[iclust][imass][imom][r],SinP[imom],SinP2[imom],r);
	  }

      

        printf("====================\n\n\n");
        

	////////////////////////////////
	// INVERTING		////////
	////////////////////////////////

        printf("====================\n");
	printf("INVERTING PROPAGATOR\n");

        colorcolorspinspin ****s0FT_inv=malloc(sizeof(colorcolorspinspin***)*(Nclust+1));
        for (int iclust=0; iclust<Nclust+1; iclust++) s0FT_inv[iclust]=malloc(sizeof(colorcolorspinspin**)*(Nmass));
        for (int iclust=0; iclust<Nclust+1; iclust++)  for(int imass=0; imass<Nmass; imass++) s0FT_inv[iclust][imass]=malloc(sizeof(colorcolorspinspin*)*(Nmom));
	for (int iclust=0; iclust<Nclust+1; iclust++)  for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++) s0FT_inv[iclust][imass][imom]=malloc(sizeof(colorcolorspinspin)*(Nr));



	complex **mdum=(complex**)malloc(sizeof(complex)*(12));
	for(int i=0; i<12; i++) mdum[i]=(complex*)malloc(sizeof(complex)*(12));

        complex **idum=(complex**)malloc(sizeof(complex)*(12));
        for(int i=0; i<12; i++) idum[i]=(complex*)malloc(sizeof(complex)*(12));

	colorcolorspinspin s0_dump;

	int ind_ext, ind_int;
	int iflag=0; 

	for (int r=0; r<Nr; r++)
        for (int imom=0; imom<Nmom; imom++)
        {
	    for (int iclust=0; iclust<Nclust+1; iclust++) for (int imass=0; imass<Nmass; imass++)
            {
                        for (int ic_ext=0; ic_ext<3; ic_ext++) for (int is_ext=0; is_ext<4; is_ext++) for (int ic_int=0; ic_int<3; ic_int++) for (int is_int=0; is_int<4; is_int++)
                         {
				ind_ext=is_ext+4*(ic_ext);
				ind_int=is_int+4*(ic_int);
				mdum[ind_int][ind_ext]=s0FT[iclust][imass][imom][r][ic_int][ic_ext][is_int][is_ext];			

			}
             		inverse(mdum,12,idum,iflag);
			if (iflag==1) {break; printf("Light quark propagator is singular\n");}
			for (int ic_ext=0; ic_ext<3; ic_ext++) for (int is_ext=0; is_ext<4; is_ext++) for (int ic_int=0; ic_int<3; ic_int++) for (int is_int=0; is_int<4; is_int++)
                         {
                                ind_ext=is_ext+4*(ic_ext);
                                ind_int=is_int+4*(ic_int);
				s0FT_inv[iclust][imass][imom][r][ic_int][ic_ext][is_int][is_ext]=idum[ind_int][ind_ext];
			}
			/////////////////////////
			// Inversion Test
			///////////////////////
			for (int ic_ext=0; ic_ext<3; ic_ext++) for (int is_ext=0; is_ext<4; is_ext++) for (int ic_int=0; ic_int<3; ic_int++) for (int is_int=0; is_int<4; is_int++)
			{
				s0_dump[ic_int][ic_ext][is_int][is_ext]=0.;
				for (int ic_dump=0; ic_dump<3; ic_dump++) for (int is_dump=0; is_dump<4; is_dump++)
						s0_dump[ic_int][ic_ext][is_int][is_ext]=s0_dump[ic_int][ic_ext][is_int][is_ext]+s0FT[iclust][imass][imom][r][ic_int][ic_dump][is_int][is_dump]*s0FT_inv[iclust][imass][imom][r][ic_dump][ic_ext][is_dump][is_ext];
			}
			for (int ic_ext=0; ic_ext<3; ic_ext++) for (int is_ext=0; is_ext<4; is_ext++) for (int ic_int=0; ic_int<3; ic_int++) for (int is_int=0; is_int<4; is_int++)
			{
				if (ic_ext != ic_int || is_ext != is_int) 
				{
					if (cabs(s0_dump[ic_int][ic_ext][is_int][is_ext])> 0.1) 
						printf("Inversion failed for imom=%d (s0*s0inv)_%d%d%d%d=%8.5f+I%8.5f \n",imom,ic_int,ic_ext,is_int,is_ext,creal(s0_dump[ic_int][ic_ext][is_int][is_ext]),cimag(s0_dump[ic_int][ic_ext][is_int][is_ext]));
				}
				else
                                 {
				       if (cabs(s0_dump[ic_int][ic_ext][is_int][is_ext])-1> 0.1) 
						printf("Inversion failed  for imom=%d (s0*s0inv)_%d%d%d%d=%8.5f+I%8.5f \n",imom,ic_int,ic_ext,is_int,is_ext,creal(s0_dump[ic_int][ic_ext][is_int][is_ext]),cimag(s0_dump[ic_int][ic_ext][is_int][is_ext]));
				}
				
			}
                        ////////////////////////////////////
                        //  end inversion test
                        ////////////////////////////////////

	     }
        }
	printf("====================\n\n\n");

       /////////////////////////////////////////////////////
        // FORM FACTOR OF THE INVERTED FOURIER-TRANSFORMED PROPAGATOR
        //////////////////////////////////////////////////////

	printf("====================\n");
        printf("FORM FACTOR OF THE INVERTED PROPAGATOR\n");
        complex *****s0_inv_ff=malloc(sizeof(complex****)*(Nclust+1));
        for (int iclust=0; iclust<Nclust+1; iclust++) s0_inv_ff[iclust]=malloc(sizeof(complex***)*(Nmass));
        for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) s0_inv_ff[iclust][imass]=malloc(sizeof(complex**)*(Nmom));
	for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++) s0_inv_ff[iclust][imass][imom]=malloc(sizeof(complex*)*(Nr));
        for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++) for (int r=0; r<Nr; r++) s0_inv_ff[iclust][imass][imom][r]=malloc(sizeof(complex)*(3));


       for (int r=0; r<Nr; r++)
       for (int imom=0; imom<Nmom; imom++)
       {
		for (int iclust=0; iclust<Nclust+1; iclust++) for (int imass=0; imass<Nmass; imass++)
                compute_s0tw_inv_ff(s0FT_inv[iclust][imass][imom][r],s0_inv_ff[iclust][imass][imom][r],SinP[imom],SinP2[imom],r);
       }



    


	printf("====================\n\n\n");
  
        /////////////////////////////////////////////////////
        // Z_q=-i/12 * Tr(sum_mu gamma_mu*SinP(mu)*S^-1(p))/SinP2
        //////////////////////////////////////////////////////

        printf("====================\n");
        printf("Z_q\n");

        int idir;
	double Pmu;

        complex ****Zq=malloc(sizeof(complex***)*(Nclust+1));
        for (int iclust=0; iclust<Nclust+1; iclust++) Zq[iclust]=malloc(sizeof(complex**)*(Nmass));
        for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) Zq[iclust][imass]=malloc(sizeof(complex*)*(Nmom));
	for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++) Zq[iclust][imass][imom]=malloc(sizeof(complex)*(Nr));

	for (int r=0; r<Nr; r++)
        for (int iclust=0; iclust<Nclust+1; iclust++) for (int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++)
	{
			Zq[iclust][imass][imom][r]=0.;
			for (int ic=0; ic<3; ic++)  for (int is=0; is<4; is++)
			{
			for (int mu=1; mu<5; mu++)
			{
                                if (mu==4) idir=0; else idir=mu;        
				Pmu=SinP[imom][idir];
				Zq[iclust][imass][imom][r]=Zq[iclust][imass][imom][r]-(I*Pmu*((base_gamma[mu]).entr[is])*s0FT_inv[iclust][imass][imom][r][ic][ic][(base_gamma[mu]).pos[is]][is])/(12*SinP2[imom]);
				
			}
			}
	}
       printf("====================\n\n\n");

       /////////////////////////////////////////////////////
       // Lambda_O'= Tr [ G_amp_O' * P_O'] ; O' = e^(-ig5 pi/4) O e(ig5 pi/4) 
       // Zq^-1 Z_O * Lambda_O' = 1 
       // The appeareance of O' is a trivial consequence of the fact that RCs,  are named after the (here unphysical) twisted quark basis while the operator OΓ = uΓd is expressed in the physical quark basis.
       // 
       // Lambda= Tr [ G_amp_O' * P_O' ] ; 
       // Twisted Mass combinations : r1=-r2 --> u,d
       // G_amp_O' = S_u^-1 G_O' S_d^-1 = S_u^-1 G_O' g5 S_u^\dagger g5
       // G_O'=<u(x) O' \bar d(x) > =S_u O' S_d^back=S_u O' g5 S_u^\dagger g5
       // Using TM fermions ZA <-> ZV
       // BUT if instead O' I use O both effects compensate: A'=-iV; V'=iA; P'=P; S'=S; T'=T
	// and i factors compensate with P_O=conj(O)

	//////////////////////////////////////////////////
	// So I use Lambda_O and I then I can read directly ZA and ZV
	/////////////////////////////////////////////////

       // P_O projector onto tree-level
       // if O=I,g5,g_mu,g5*g_mu... P_O=conj(O)
       //////////////////////////////////////////////////////

        printf("====================\n");
        printf("LAMBDA AND Z_O\n");

      complex *****Lambda=malloc(sizeof(complex****)*(16));
      for (int igamma=0; igamma<16 ; igamma++) Lambda[igamma]=malloc(sizeof(complex***)*(Nclust+1));
      for (int igamma=0; igamma<16 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) Lambda[igamma][iclust]=malloc(sizeof(complex**)*(Nmass*Nmass));
      for (int igamma=0; igamma<16 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass*Nmass; imass++) Lambda[igamma][iclust][imass]=malloc(sizeof(complex*)*(Nmom));
      for (int igamma=0; igamma<16 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass*Nmass; imass++) for (int imom=0; imom<Nmom; imom++)  Lambda[igamma][iclust][imass][imom]=malloc(sizeof(complex)*(Nr));

      complex *****ZO=malloc(sizeof(complex****)*(16));
      for (int igamma=0; igamma<16 ; igamma++) ZO[igamma]=malloc(sizeof(complex***)*(Nclust+1));
      for (int igamma=0; igamma<16 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) ZO[igamma][iclust]=malloc(sizeof(complex**)*(Nmass*Nmass));
      for (int igamma=0; igamma<16 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass*Nmass; imass++) ZO[igamma][iclust][imass]=malloc(sizeof(complex*)*(Nmom));
      for (int igamma=0; igamma<16 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass*Nmass; imass++) for (int imom=0; imom<Nmom; imom++) ZO[igamma][iclust][imass][imom]=malloc(sizeof(complex)*(Nr));

      // auxiliares (solo por claridad de las formulas)

      colorcolorspinspin *G_r=malloc(sizeof(colorcolorspinspin)*(Nconf));
      colorcolorspinspin *G=malloc(sizeof(colorcolorspinspin)*(Nclust+1));
      colorcolorspinspin S_xy, S_yx, Sinv_xy, Sinv_yx;
      int imass;

      for (int igamma=0; igamma<16; igamma++)
      {
      whichGamma(igamma,gammaprint);
      printf("%s\n",gammaprint);
      for (int r=0; r<Nr; r++)
      for (int imom=0; imom<Nmom; imom++)
      {
	   for (int im2=0; im2<Nmass; im2++) for (int im1=0; im1<Nmass; im1++)
	   {
		imass=im2*Nmass+im1;
		//////////////////////////////////////////////
                // Compute NON amputed Green Function
		/////////////////////////////////////////////

		// On one single configuration the quark propagator is not transationally invariant. It is only recovered after averaging over the configurations
		// Then we compute G over the initial configurations 
		// and then we average
                // G(p)=\sum_n G_n(p)/N
		for (int iconf=0; iconf<Nconf; iconf++) 
	        {
			// 1. Initialize
			for (int ic1=0; ic1<3 ; ic1++) for (int ic2=0; ic2<3 ; ic2++) 
                        for (int is1=0; is1<4 ; is1++) for (int is2=0; is2<4 ; is2++) 
			{
				S_xy[ic1][ic2][is1][is2]=s0FT_r[iconf][im1][imom][r][ic1][ic2][is1][is2];
				S_yx[ic1][ic2][is1][is2]=((base_gamma[5]).entr[is1])*conj(s0FT_r[iconf][im2][imom][r][ic2][ic1][(base_gamma[5]).pos[is2]][(base_gamma[5]).pos[is1]])*(base_gamma[5]).entr[is2];
			}
			// 2. Compute
                       for (int ic1=0; ic1<3 ; ic1++) for (int ic2=0; ic2<3 ; ic2++)
                       for (int is1=0; is1<4 ; is1++) for (int is2=0; is2<4 ; is2++)
			{
				G_r[iconf][ic1][ic2][is1][is2]=0.;
				for (int ic3=0; ic3<3 ; ic3++) for (int is3=0; is3<4 ; is3++)
					G_r[iconf][ic1][ic2][is1][is2]=G_r[iconf][ic1][ic2][is1][is2]+S_xy[ic1][ic3][is1][is3]*((base_gamma[igamma]).entr[is3])*S_yx[ic3][ic2][(base_gamma[igamma]).pos[is3]][is2];
			}
		}
	
		//3. Clusterize
                for (int ic1=0; ic1<3 ; ic1++) for (int ic2=0; ic2<3 ; ic2++)
                for (int is1=0; is1<4 ; is1++) for (int is2=0; is2<4 ; is2++)
                {
                       for(int iconf=0;iconf<Nconf;iconf++) aux_conf[iconf]=G_r[iconf][ic1][ic2][is1][is2];
                       conf2clust(aux_conf, aux_clust,Nconf,Nclust,Nex);
                       for(int iclust=0;iclust<Nclust+1;iclust++) G[iclust][ic1][ic2][is1][is2]=aux_clust[iclust];
                }

		/////////////////////////////
                // Compute Lambda and Z_O
		///////////////////////////
		 for (int iclust=0; iclust<Nclust+1; iclust++)
                {

			// 1. Initialize
			Lambda[igamma][iclust][imass][imom][r]=0.;
			for (int ic1=0; ic1<3 ; ic1++) for (int ic2=0; ic2<3 ; ic2++)
                        for (int is1=0; is1<4 ; is1++) for (int is2=0; is2<4 ; is2++)
                        {
				Sinv_xy[ic1][ic2][is1][is2]=s0FT_inv[iclust][im1][imom][r][ic1][ic2][is1][is2];
				Sinv_yx[ic1][ic2][is1][is2]=((base_gamma[5]).entr[is1])*conj(s0FT_inv[iclust][im2][imom][r][ic2][ic1][(base_gamma[5]).pos[is2]][(base_gamma[5]).pos[is1]])*(base_gamma[5]).entr[is2];
			}
			// 2. Compute
			for (int ic1=0; ic1<3 ; ic1++) for (int ic2=0; ic2<3 ; ic2++) for (int ic3=0; ic3<3 ; ic3++)
			for (int is1=0; is1<4 ; is1++) for (int is2=0; is2<4 ; is2++) for (int is3=0; is3<4 ; is3++)
			{
			Lambda[igamma][iclust][imass][imom][r]=Lambda[igamma][iclust][imass][imom][r]+Sinv_xy[ic1][ic2][is1][is2]*G[iclust][ic2][ic3][is2][is3]*Sinv_yx[ic3][ic1][is3][(base_gamma[igamma]).pos[is1]]*conj((base_gamma[igamma]).entr[is1])/12.;
			}

			// Z0=Zq/Lambda
			ZO[igamma][iclust][imass][imom][r]=sqrt(Zq[iclust][im1][imom][r]*Zq[iclust][im2][imom][r])/Lambda[igamma][iclust][imass][imom][r];

		}
	
	    }	

      }
     }

     printf("====================\n\n\n");





///***************************************************************************************************************///////
//********************************************************************************************************************///
///***************************************************************************************************************///////

	/////////////////
	// Write output
	////////////////



     complex **aux_allmomenta=malloc(sizeof(complex*)*(Nmom));
     for (int imom=0; imom<Nmom; imom++) aux_allmomenta[imom]=malloc(sizeof(complex)*(Nclust+1));

      
      printf(" WRITING OUTPUT IN %s\n",dir_output);
      for (int sigma=0; sigma<3; sigma++)
      {

        for (int r=0; r<Nr; r++) for (int imass=0; imass<Nmass; imass++)
        {
        // ff S(p)

         sprintf(s0ff_output,"%s_r%d_m%d_sigma%d.out",s0ff_base_output,r,imass,sigma+1);
         for (int imom=0; imom<Nmom; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) aux_allmomenta[imom][iclust]=s0_ff[iclust][imass][imom][r][sigma];
         write_output_allmomenta(s0ff_output,aux_allmomenta,SinP2);

         sprintf(file,"ff_r%d_m%d_sigma%d",r,imass,sigma+1);
         write_output_momentabymomenta(dir_output,file,aux_allmomenta,SinP2);


        }
      }
      for (int sigma=0; sigma<3; sigma++)
      {

        for (int r=0; r<Nr; r++) for (int imass=0; imass<Nmass; imass++)
        {

        // ff S^-1(p)
         sprintf(s0ff_output,"%s_inv_r%d_m%d_sigma%d.out",s0ff_base_output,r,imass,sigma+1);
         for (int imom=0; imom<Nmom; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) aux_allmomenta[imom][iclust]=s0_inv_ff[iclust][imass][imom][r][sigma];
         write_output_allmomenta(s0ff_output,aux_allmomenta,SinP2);

         sprintf(file,"ff_inv_r%d_m%d_sigma%d",r,imass,sigma+1);
         write_output_momentabymomenta(dir_output,file,aux_allmomenta,SinP2);

        }
      }

    //////////////
    // Zq
    ////////////////////
    for (int r=0; r<Nr; r++) for (int imass=0; imass<Nmass; imass++)
     {
	 sprintf(z_output,"%sq_r%d_m%d.out",z_base_output,r,imass);
         for (int imom=0; imom<Nmom; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) aux_allmomenta[imom][iclust]=Zq[iclust][imass][imom][r];
         write_output_allmomenta(z_output,aux_allmomenta,SinP2);

         sprintf(file,"Zq_r%d_m%d",r,imass); 
         write_output_momentabymomenta(dir_output,file,aux_allmomenta,SinP2);
     }


     /////////////////////
    // ZO 
    //////////////////// 
   for (int igamma=0; igamma<16; igamma++)
   {
	whichGamma(igamma,gammaprint);
	for (int r=0; r<Nr; r++)
        for (int im2=0; im2<Nmass; im2++) for (int im1=0; im1<Nmass; im1++)
         { 
         imass=im2*Nmass+im1;

         sprintf(z_output,"%s%s_r%d_m%d%d.out",z_base_output,gammaprint,r,im1,im2);
         for (int imom=0; imom<Nmom; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) aux_allmomenta[imom][iclust]=ZO[igamma][iclust][imass][imom][r];
         write_output_allmomenta(z_output,aux_allmomenta,SinP2);

         sprintf(file,"Z%s_r%d_m%d%d",gammaprint,r,im1,im2);                
         write_output_momentabymomenta(dir_output,file,aux_allmomenta,SinP2);


       }
    }

     /////////////////////
    // Lambda 
    //////////////////// 
   for (int igamma=0; igamma<16; igamma++)
   {
        whichGamma(igamma,gammaprint);
	for (int r=0; r<Nr; r++)
        for (int im2=0; im2<Nmass; im2++) for (int im1=0; im1<Nmass; im1++)
        {
          imass=im2*Nmass+im1;

         sprintf(Lambda_output,"%s%s_r%d_m%d%d.out",Lambda_base_output,gammaprint,r,im1,im2);
         for (int imom=0; imom<Nmom; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) aux_allmomenta[imom][iclust]=Lambda[igamma][iclust][imass][imom][r];
         write_output_allmomenta(Lambda_output,aux_allmomenta,SinP2);

         sprintf(file,"Lambda%s_r%d_m%d%d",gammaprint,r,im1,im2);
         write_output_momentabymomenta(dir_output,file,aux_allmomenta,SinP2);


       }
    }
     



   //all Lambda together
   for (int r=0; r<Nr; r++)
   for (int im2=0; im2<Nmass; im2++) for (int im1=0; im1<Nmass; im1++)
   {
       imass=im2*Nmass+im1;

	for (int imom=0; imom<Nmom; imom++)
	{

             sprintf(Lambda_output,"%s/p%03d/Lambda_r%d_m%d%d.dat",dir_output,imom,r,im1,im2);
             FILE *fout=fopen(Lambda_output,"w");
             if(fout == (FILE *)NULL) {
                        fprintf(stderr,"Unable to open file %s for writing\n",Lambda_output);
                        return EXIT_FAILURE;
             }
	     for (int igamma=0; igamma<16; igamma++)
	     {
                 for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]=Lambda[igamma][iclust][imass][imom][r];
                 fprintf(fout,"%8.10f %8.10f %8.10f \n",SinP2[imom],creal(Lambda[igamma][Nclust][imass][imom][r]),creal(error_jack(aux_clust,Nclust)));
	     }
             fclose(fout);
        }
   }



     printf("====================================================\n\n\n");

////// FINISH
	printf("Free memory\n");

	free(s0_ff); free(s0_inv_ff);
	free(Zq); free(Lambda); free(ZO);
	
	printf("Calc_FF_AmpGF_and_Z.c finished OK ! \n");
	
  return 0;
}
